Evidence of a spin liquid phase in the frustrated honeycomb lattice 
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In the present paper we present some new data supporting the existence of a spin-disordered phase 
in the Heisenberg model on the honeycomb lattice with antiferromagnetic interactions up to third 
neighbors along the line J2 = J3, predicted in^. We use the Schwinger boson technique followed 
by a mean field decoupling and exact diagonalization for small systems to show the existence of an 
intermediate phase with a spin gap and short range Neel correlations in the strong quantum limit 
(S=i). 



PACS numbers: 



o 



I 

o 
o 



> 
00 
00 
00 
vn 

cn 
o 



INTRODUCTION 



Geometrical frustration in two-dimensional (2D) anti- 
ferromagnets is expected to enhance the effect of quan- 
tum spin fluctuations and hence suppress magnetic order 
giving rise to a spin liquid^ and this idea has motivated 
many researchers to look for its realization^ri. 

One candidate to test these ideas is the honeycomb 
lattice, which is bipartite and has a classical Neel ground 
state, but due to the small coordination number [z — Z), 
quantum fluctuations could be expected to be stronger 
than those in the square lattice and may destroy the an- 
tiferromagnetic long-range order (LRO}^i^. 

The study of frustrated quantum magnets on the hon- 
eycomb lattice has also experimental motivations"'^^"^'*. 

The analysis of the hexagonal lattice from a more gen- 
eral point of view has gained lately a lot of interest both 
coming from graphene-related issues and from the possi- 
ble spin-liquid phase found in the Hubbard model in such 
geometrjJ^""— . 

Motivated by the previous results, in this paper we 
show the study of the Heisenberg model on the honey- 
comb lattice with first (Ji), second (J2) and third (J3) 
neighbors coupling&i§,^ along the special line J2 = J3. 
Using Schwinger boson mean field theory (SBMFT) and 
exact diagonalization we find strong evidence for the ex- 
istence of an intermediate disordered region where a spin 
gap opens and spin-spin correlations decay exponentially. 
Using exact diagonalization of small clusters we also have 
calculated the dimer-dimer correlation function that in- 
dicates short range dimer-dimer order. Although our re- 
sults correspond to a specific line, we conjecture that 
the quantum disordered phase that we have found in the 
vicinity of the tricritical point extends within a finite re- 
gion around it. Previous evidence of massive behaviour in 
the hexagonal lattice Heisenberg model has been found 
in other regions of the phase space by means of exact 
diagonalization in ref J^. 

The Heisenberg model on the Ji — J2 — J3 honeycomb 




FIG. 1: (Color online) Honeycomb lattice with sublattices A 
and B. The vectors ei = (-^, |) and 6*2 = (-^, -|) are the 
primitive traslation vectors of the direct lattice. 



lattice is given by 

= Jl ^ Sx ■ Sy + J2 ^ Sx ■ Sy + J3 ^ Sx ■ Sy, (l) 

<xy)i <''y>2 (xy>3 

where Sx is the spin operator on site x and (xy)„ in- 
dicates sum over the n-th neighbors (see Fig. [1]). In 
the classical limit, S — t- 00, the model displays differ- 
ent zero temperature phases with a tricritical point at 
J2 — J3 — ^Ji- At this particular point the classical 
ground state has a large GS degenerac y-^^'^^ . The Heisen- 
berg model on the honeycomb lattice was studied using 
SBMFT by Mattsson et al^ for antiferromagnetic inter- 
actions at first and second neighbors. Here we study the 
Hamiltonian ([T} using a rotationally invariant version of 
this technique, which has proven successful in incorpo- 
rating quantum fluctuations^Sr— . 



2. SCHWINGER BOSONS MEAN-FIELD 
THEORY. 

In this section we describe in detail the Schwinger bo- 
son mean field theory used in the present work. The 
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SU{2) Heisenberg Hamiltonian on a general lattice can 
be written as 



(2) 



xya/3 



where x and y are the positions of the unit cells and vec- 
tors Ta correspond to the positions of each atom within 
the unit cell. Jap{x — y) is the exchange interaction be- 
tween the spins located in x -|- e y -|- r^. 

In what follows we assume that the classical order can 
be parameterized as 



QX 

'-'x+ro 

qy 

•^x+rc 



= S'siny>a(x) 
= 

= S cos ipa{x), 



(3) 
(4) 
(5) 



with i/5„ (x) = Q • X + where Q is the ordering vector 
and 9a are the relative angles between the classical spins 
inside each unit cell. 

The spin operators Sx on site x are represented by two 
bosons 6x<7 =t)4-) 



Sr=-bt.a-br 



(6) 



where a = (cj;, CTj,, cr^) arc the Pauli matrices. This is 
a faithful representation of the algebra SU(2) if we take 
into account the following local constraint 



2S = 6l,Sxt 



Km- 



(7) 



In this representation, the exchange term can be ex- 
pressed as 



where Aa.j3{x, y) and Ba,p{^, y) are the SU{2) invariants 
defined as 



(T 

B„,^(x,y) = \ Y.bi(^^% (10) 



with a =t.i and double dots (: O :) indicate normal or- 
dering of operator O. This decoupling is particularly use- 
ful to the description of magnetic systems near disordered 
states, because it allows to treat antiferromagnetism and 
ferromagnetism in equal footing. 

To construct a mean field theory we perform a Hartree- 
Fock decoupling 



(S, 



with 



^a/3(x- y) 
-B;^/3(x-y) 

{{^x+ra^y+r b)mf) 



[S*;3(x-y)Bc,;3(x, y) 

^a/3(x - y)i„,3(x, y) -f H.c] 

{(Sx+rc ■ Sy+rajMF), (H) 



(il,(x,y)> (12) 
(Sl^(x,y)) (13) 
\B^p{x-y)\^ -\A^p{x-y)\\ 



Sx-l-r 



Sy+r^ = : -B^«(x,y)Bc.,3(x,y) : -iL(x,y)i„^ 



and where ( ) denotes the expectation value in the ground 
state at T = 0. Because several functions involved in the 
Hamiltonian depend on the difference x — y we change 
^aJriables to R = x — y and eliminating x in the sums we 
()$))tain 



Hmf = \Y. ■^"'^(R) j i E K.MR) b'^+y.M - <yA^A^) bli?y,.6;^,a + H.c] - ( \B^Am'' - \Ao.Am') \ • (14) 



The mean field Hamiltonian is quadratic in the boson 
operators and can be diagonalized in real space. How- 
ever, as we look for translational invariant solutions, it 
is convenient to transform the operators to momentum 
space 



S(«) 



1 



(a) ik;-(x+r„) 



(15) 



After some algebra and using the symmetry properties: 

(16) 

we obtain the following form for the Hamiltonian 



Ja0(R) 
A„/3(R) 

-Bo/3(R) 



J/3a(-R) 

-A^„(-R) 
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k:a/3 cr 

E w [i^-z'Wi' - i^^^wfl . 

Ra/3 



(17) 



where 



B 



A 



(k) = \ ^J„^(R)i3„;3(R)e-*-(^+--^) (18) 

R 

» = i;^J.^(R)A<,;3(R)e- 



7^/3(k) = i^Jc^(R)A<,^(R)c- 



ik.{R+r„-r^) 



ik.{R+r„-rg) 



(19) 
(20) 



with 



xa \ <y / 



(22) 



Now, we impose the constraint ([7]) in average over each 
sublattice a by means of Lagrange multipUers A*-"-* 



Hmf ^ Hmf + H\ 



Using the symmetries ([T6|) we can see that both kinds 
of bosons (t, i) give the same contribution to the Hamil- 
(21) Ionian. Then, we can perform the sum over a to obtain 



J 



kQ/3 



Rq/3 



It is convenient to define a vector operator b^^(k) — 
^bj^^,b-k;) where 



'kt ~ v'^kt ' "kf 



^kt 



(23) 

b-ki = {Vit^:h%{,.-.:b'^tf) (24) 

and ric is the number of atoms in the unit cell. Now, we 
can write the Hamiltonian as 

Hmf = E bt(k).i?(k).b(k) (25) 

k 

- (2^+l)iV,E^'"^ - {Hmf), 



Z?(k) 



where the 2nc x 2nc dynamical matrix 13 (k) is given by 

'7f^(k) + A("'5„^ -7^^(k) 

7^;3(k) 7f^(k) + A('^'5„^ 

To diagonalize the Hamiltonian (|25p we need to per- 
form a para-unitary transformation of the matrix -D(k) 
that preserves the bosonic commutation relations. We 
can diagonalize the Hamiltonian by defining the new op- 
erators a = F • b, where the matrix F satisfy 



(Ft)-i . ra • {F)-^ = rs. 



T3 



^2x2 
-l2> 



(26) 



With this transformation, the Hamiltonian reads 
Hmf = E 4 • E(k) • at - 2{S + l)iV, E A^"' - {H 



MF 



(27) 



where 



E(k) =diag(wi(k),...,cj„(k),cji(k),...,w„(k)). (28) 

In term of the original bosonic operators, the mean 
field parameters are 

k 
k 

- e-'''^^+''--'^^\b^}[lb%)} (30) 

and the constraint in the number of bosons can be written 
in the momentum space as 



= 2SNc 



(31) 
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FIG. 2: (a): GS energy per unit cell for S = | as a function 
of J2/J1 for a lattice of 32 sites. The circles are exact results 
(ED) and the squares are the SBMFT results, (b): Spin- 
Spin correlation function (SSCF) vs distance X in the zig-zag 
direction obtained within SBMFT for a 32 x 32 system 
For < J2/J1 < 0.41, the SSCF correspond to the Neel 
phase with long-rage-order (LRO), for 0.41 < J2/J1 < 0.6 
the correlations are short ranged indicating a gap zone with 
sort-range-order (SRO), and for 0.6 < J2/J1 the correlations 
correspond to the collinear phase (ferromagnetic correlations 
in the zig-zag direction). The inset in Fig. (b) shows the 
finite size results for the SSCF obtained by ED and SBMFT 
for 32 sites. Lines are guides to the eye and different colors 
are used for clarity. 



where Nc is the total number of unit cells and S is the 
spin strength. The mean field equations ((29)) and (|30)) 
must be solved in a self-consistent way together with the 
constraints pil) on the number of bosons. Finding nu- 
merical solutions involves finding the roots of 24 coupled 
nonlinear equations for the parameters A and B, plus the 
additional constraints to determine the values of the La- 
grange multipliers A*^"-*. We perform the calculations for 
finite but very large lattices and finally we extrapolate 
the results to the thermodynamic limit. We solve nu- 
merically for several values of the frustration parameter 
J2/J1 and with the values obtained for the MF parame- 
ters and the Lagrange multipliers we compute the energy 
and the new values for the MF parameters. We repeat 
this self-consistent procedure until the energy and the 
MF parameters converge. After reaching convergence we 
can compute all physical quantities like the energy, the 
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FIG. 3: Gap in the boson dispersion as a function of J2/J1 
for S = 1/2 from Ref.— . In the region J2/J1 ~ 0.6 the system 
remains gapped. Inset: finite size scaling for the gap. (i) 
J2/J1 — 0.5 (7 — 0.6451), (ii): Circles correspond to J2/J1 = 
0.05 (7 = 0.911) and squares correspond to J2/J1 = 0.35 
(7 = 0.758). 



spin-spin correlations and the excitation gap. In order 
to support the analytical results of the MF approach, we 
have also performed exact diagonalization on finite sys- 
tems with 18, 24 and 32 spins with periodic boundary 
conditions for S = 1/2 using Spinpack^"^. 



3. RESULTS 

In Fig. HJa) we show the ground state energy per unit 
cell as a function of the frustration for a system of 32 
sites calculated by means of SBMFT and ED, showing 
an excellent agreement between both approaches. The 
advantage of the SBMFT is that it allows to study much 
larger systems: we have studied different system sizes 
up to 3200 sites and extrapolated to the thermodynamic 
limit. 

For the present model we only find commensurate 
collinear phases and for these phases, the wave vector 
Qo = Q/2 where the dispersion relation has a minimum 
remains pinned at a commensurate point in the Brillouin 
zone, independently of the value of the frustration J2/J1. 
In the thermodynamic limit, a state with LRO is char- 
acterized in the Schwinger boson approach by a conden- 
sation of bosons at the wave vector Qq. This implies 
that the dispersion of the bosons in a state with LRO is 
gapless. As we discussed earlier, we solve the mean field 
equations for finite systems, then to detect LRO we cal- 
culate the gap in the bosonic spectrum as a function of 
J2/ Ji for different system sizes and perform a finite size 
scaling finding a finite region where the system remains 
gapped. 

The structure of the different phases can be understood 
calculating the spin-spin correlation function (SSCF). 
For J2/J1 < 0.41 the SSCF is antiferromagnetic in all 
directions and is long-ranged while for 0.6 < J2/J1 
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FIG. 4: (Color online) Inverse of the critical spin Sc as a 
function of J2/J1 obtained using SBMFT by us in^. For the 
5=1/2 case, there is a range 0.41 < J2/J1 < 0.6 where the 
system has a spin-gap indicating a quantum disordered phase 
(see Fig. [3|. The dotted- line correspond to the classical limit 
S — >■ 00 where the ground state correspond to the Neel phase 
with Q = (0, 0) and (/>a — i^s = tt in the region J2/J1 < 0.5, 
while for J2/J1 > 0.5 the ground state correspond to the CAF 
phase characterized by Q = (27r/-\/3, 0) and (j)A — (I>b ~ i"- 
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FIG. 5: Dimer-dimer correlations D{i,2),(fc,i) between the ref- 
erence bond (1, 2) and bonds (k, I) in the ground-state of the 
iV = 24 sample for J2/J1 ~ 0.55. The number on bond {k,l) 



indicates the value of Dn 



{l,2),{fc,0 



truncated to the two first 



significant digits. Full black (orange) lines indicate positive 
(negative) values of I5(i,2),(fe,o • 



we have found ferromagnetic LRO correlations in the 
zig-zag direction that correspond to the CAF phase. 
The most interesting result is in the intermediate re- 
gion 0.41 < J2/J1 < 0.6 where the resuhs for the SSCF 
predict a quantum disordered state with a gap in the 
bosonic dispersion and the spin-spin correlation function 
shows Neel short range order. A plot of the SSCF for 
J^/J^ = 0,0.2,0.5 and 0.8 obtained within SBMFT is 
presented in Fig. [^b- In Fig. S] we show the ground state 



phase diagram as a function of The classical phase 

diagram reduces to that shown in the line 1/S'c = of 
Fig. H] where two collinear phases meet at the classical 
critical point J2/J1 = 0.5. On the one hand, for 1/S 
smaller than a critical value l/Sc{J2/Ji), the correlation 
functions have LRO, characterized by a condensation of 
bosons at the wave vector Qq. On the other hand, when 
1/S is greater than l/Sc{J2/ Ji), the correlation func- 
tions have SRO indicating quantum disorder. 
In the intermediate region the results found with SBMFT 
predict a quantum disordered region 0.41 < J2/J1 < 0.6. 
In this region a gap opens in the bosonic dispersion and 
the spin-spin correlation function shows Neel short range 
order followed by the LRO CAF phase for J2/J1 > 0.6. 
In Fig. [3] we show the extrapolation of the boson gap 
as a function of the frustration. The inset shows an ex- 
ample of the finite size scaling for different values of the 
frustration. 

Previous results show that for 0.41 < J2/J1 < 0.6 the 
ground state has no magnetic order—. The main question 
now is: Is this non magnetic quantum phase a quantum 
disordered one? Or does it exhibit any other kind of non- 
magnetic order?. To answer this question the knowledge 
of the spin-spin correlation function is not enough and 
one has to check for other types of (non-magnetic) order- 
ing patterns. 

One kind of non magnetic order that can set in is the 
dimer long-range order. The dimer operator on a pair of 
sites is defined as dij = 1/4 — • Sj, and one usu- 
ally defines the dimer-dimer correlation between bonds 
(ij) and (kj) as L»(jj)_(fc_j) =< dt,jdk,i > - < dij >< 
dfe_/ >. In order to understand the nature of the ground 
state in the intermediate region, we have calculated de 
dimer-dimer correlation function defined above by means 
of exact diagonalization on a 24 sites cluster with peri- 
odic boundary conditions for S — 1/2. The correlation 
pattern for dimers on first neighbor bonds is displayed in 
Fig. [5l We can see that the exact dimer-dimer correla- 
tions show a rather fast decay suggesting that there is no 
dimer order in the groud state, though due to the small 
size of the cluster studied, this is not conclusive and we 
cannot discard other ordering patterns. 

In summary, the results presented here further support 
the existence of a region in the intermediate frustration 
regime where the system does not show quantum mag- 
netic order for S = 1/2. 

Note added: When this manuscript was completed 
we became aware of two independent works providing 
an analysis of the model using a combination of ex- 
act diagonalizations^ii^ and an effective quantum dimer 
model, as well as a self-consistent cluster mean-field 
theorj^. Several similar findings show a good corre- 
spondence of both approaches. 
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